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Summary: 

Research  was  conducted  to  evaluate  the  feasibility  of  wavelet  transforming 
second  order  statistics.  To  identify  the  modulation  type  and  to  extract  signal 
modulation  parameters  of  frequency  hopped  signals  the  wavelet  transform  is 
applied  to  the  2-dimensional  instantaneous  correlation  function.  Parameters  of 
interest  are  switching  times  and  hop  frequencies.  The  characteristics  of  the 
wavelet  transform  of  the  correlation  function  are  derived.  The  wavelet  transform 
of  the  correlation  function  can  replace  the  requirement  to  Fourier  transform  the 
complete  correlation  surface  with  a  Fourier  transform  once  per  hop  to  estimate  the 
hop  frequency  of  a  given  hop.  Wavelet  processing  is  applied  along  the  delay  and 
time  axes. 

Processing  along  the  delay  axis  leads  to  the  following  conclusions: 

•  To  perform  visual  identification  an  SNR  of  3  dB  is  required. 

•  For  frequencies  larger  than  1/16  of  the  sampling  rate  the  hop  frequencies  can  be 
estimated  with  a  success  rate  of  1 00%  for  SNR  levels  of  0  dB  or  better. 

•  To  estimate  the  hop  times  with  an  accuracy  of  12  to  17.5  %  an  SNR  of  6  dB  or 
better  is  required. 

•  If  the  true  start  and  stop  time  are  used  we  can  obtain  an  improvement  of  the 
spectral  estimation  performance  by  about  2  dB. 

•  Processing  along  the  time  axis  allows  detection  of  the  transition  times  of 
frequency  and  time  hopped  signals.  The  current  implementation  is  limited  to  work 
with  one  transition  per  observation  interval  but  permits  robust  detection  at  an  SNR 
level  of  3  dB. 


iii 


1.  Introduction 


This  work  focuses  on  the  wavelet  transform  of  the  second  order  moment  function  to 
enhance  detection  and  classification  of  frequency  hopped  signals.  There  are  four  useful  and 
important  representations  when  dealing  with  non-stationary  processes: 

i)  the  temporal  correlation  function, 

ii)  the  ambiguity  function, 

iii)  the  spectral  correlation  function,  and 

iv)  the  time  frequency  distribution. 

Any  of  flie  four  representations  can  be  reached  from  any  other  representation  by  performing  a 
one  or  two-dimensional  Fourier  transforms  (Fig.  1). 

We  will  use  the  temporal  correlation  function  as  the  starting  point  and  examine  the  two 
domains  reachable  by  a  one-dimensional  Fourier  transform.  If  we  Fourier  transform  over  the 
delay  or  time  axis  of  the  correlation  function,  we  obtain  the  Wigner-Ville  Distribution  or  the 
ambiguity  function,  respectively.  In  all  discussions  we  use  an  approximation  for  the  true 
correlation  function.  If  the  expectation  operator  is  left  out  then 
Rx  (tiT)  =  E[x*(t-T/2)  x(t-i-x/2)]  becomes  x*(t-T/2)  x(t-hr/2). 

This  is  not  an  imusual  r^proximation  since  we  caimot  apply  an  expectation  operator  to 
the  data  nor  can  we  apply  time  domain  averaging.  Our  work  will  focus  on  this  correlation 
function  which  is  called  the  instantaneous  correlation  function.  In  particular  we  will  rqplace  the 
traditional  Fourier  transform  with  the  wavelet  transform.  Chapters  1-2  deal  with  general 
background,  chapters  3-5  address  transformation  over  the  delay  axis,  while  chapter  6  focuses  on 
the  transformation  over  Ihie  time  axis. 


1 


2.  General  Background 

In  some  applications  there  is  the  desire  to  intercept  digital  communication  signals.  The 
task  of  intercepting  a  communication  signal  can  be  summarized  by  i)  detect  the  signal's 
presence,  ii)  classify  the  modulation  type,  iii)  estimate  the  reception  control  parameters,  iv) 
decode  the  data,  and  v)  decrypt  the  information  content.  The  process  can  be  stopped  at  any 
intermediate  step. 

Spread  spectrum  (SS)  modulation  is  a  widely  used  modulation  technique.  Frequency 
hopping  (FH)  is  a  modulation  subset  of  SS,  and  is  primarily  addressed  in  this  report. 

Many  signal  processing  tools  are  available  to  help  to  achieve  the  tasks  listed  above.  In 
particular,  correlation  processing  and  wavelet  analysis  of  the  time  domain  data  have  been  used 
for  the  interception  task.  In  this  work  we  will  address  the  merging  of  wavelet  and  correlation 
concepts  to  enhance  detection,  classification  and  signal  parameter  estimation.  For  the  interested 
reader  an  extensive  reference  section  [1-64]  is  provided. 

The  FH  signal  is  a  non-stationary  process  having  a  two-dimensional  correlation 
fimction.  Application  of  wavelet  analysis  to  correlation  fimctions  is  a  new  area  and  is  still  in  the 
exploratory  stage.  This  work  assesses  wavelet  processing  of  the  correlation  function  along  the 
delay  and  time  axis. 
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2.1  Em:RCEPTION  OF  DIGITAL  COMMUNICATION  SIGNALS 

Interception  of  communication  signals  is  of  interest  to  a  wide  range  of  applications  in 
surveillance,  intelligence,  reconnaissance,  geo-location,  spectral  monitoring  and  jamming  [1]. 
Digital  communication  systems  can  use  a  large  number  of  modulation  techniques  (i.e.,  ASK, 
BPSK,  BFSK,  QAM,  MPSK  ,MFSK ,  Spread  Spectrum).  Interception  of  digital  communication 
signals  consists  of  detection,  classification,  parameter  estimation,  decoding,  and  decryption. 

A  large  number  of  publications  address  the  interception  of  digital  communication  signals. 
Signal  processing  is  used  for  the  interception  and  can  be  grouped  into  the  following  approaches: 

•  Second  order  moments:  Spectral  analysis  and  correlation  analysis 

•  Linear:  Linear  transforms  including  the  wavelet  transform 

•  Nonlinear:  Higher-order  spectra,  spectral  correlation,  and  cyclic-feature  processing 

•  Other:  Eigen-analysis,  singular-value  decomposition,  and  stochastic  resonance. 
Demodulation  of  ASK,  QAM,  BPSK,  and  DPSK  are  addressed  in  [2-14,17,32].  Time 

domain  correlation  and  Spectral  correlation  are  used  on  FSK  related  problems  in  [2,33], 
higher-order  moments  are  used  in  [3],  and  wavelet  analysis  is  used  in  [5,17,56,65]. 

2.2  Spread  Spectrum  Signals 

Interception  of  spread  spectrum  signals  is  addressed  in  [17-31].  These  techniques  differ 
mainly  in  the  bandwidth  of  the  interception  filter(s)  relative  to  the  bandwidth  of  the  FH  signal 
and  the  number  of  parallel  channels  relative  to  the  number  of  hopping  frequencies. 

2.3  FOURIER  ANALYSIS,  TIME  FREQUENCY  DISTRIBUTIONS  AND  WAVELETS 

Signal  analysis  treats  time  signals  as  a  linear  combination  of  elementary  basis  fimctions. 
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Well-known  exsmples  are  the  Shannon,  the  Karhunen-Loeve,  the  Gram-Schniidt  expansion,  the 
Eigen-decomposition,  and  the  Fourier  analysis.  We  will  review  the  Fourier  analysis.  Time 
Frequency  Distributions  and  Wavelet  analysis  [41,52-58,61]. 

2.4  Fourier  Transform 

The  Fourier  transform  (FT)  is  the  most  popular  signal  decomposition  [36].  It  is  used  to 
decompose  stationary  signals  into  sinusoidal  or  complex  exponential  components.  A 
non-periodic  continuous  time  signal,  x(t),  can  be  represented  as 

x(t)  =  r  X(/)  df 

J  —00 

with  (1) 

=  r  x(t)  e  dt  ; 

J  — oo 

where  the  signal  and  its  transform  are  continuous  functions  of  time  and  frequency,  respectively. 

2.5  Short  Time  Fourier  Transform 

To  track  time  evolution  the  short  time  Fourier  transform  (STFT)  was  developed.  The  STFT 
windows  the  signal  around  a  given  time  instant,  performs  the  frequency  domain  analysis,  and 
repeats  the  process  at  other  time  instants  of  interest.  The  basic  assumption  is  that  the  windowed 
signal  has  a  non-time-varying  spectrum  (local  stationarity)  within  the  time  window.  The  STFT 
for  a  continuous  signal  x(t)  is  given  by 

W.'f)  =  j°°  x(t)  w  ’(/-T)  e  dt.  (2) 

where  w(t)  is  the  window  function.  X(Ct)  is  the  spectral  description  of  x(t),  ♦  denotes 
conjugation  and  the  time  window  is  centered  at  x.  If  the  window  has  a  Gaussian  sh^e,  the 
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STFT  is  called  the  Gabor  transform.  Time  and  frequency  localization  are  controlled  by  the 
effective  window  duration.  The  STFT  has  a  fixed-time  and  fixed-frequency  resolution,  which 
results  in  a  uniform  tiling  of  the  time-frequency  plane.  For  discrete  time  signals  the  STFT  is 
defined  as 


X(k,m)  =  *(”)  ® 

■where  k  = 


(3) 


2.6  Discrete  Fourier  Transform 

A  finite-length  non-periodic  discrete  signal  has  a  continuous-frequency-domain 
representation.  The  discrete  Fourier  transform  (DFT)  and  its  fast  implementation,  the  fast 
Fourier  transform  (FFT),  use  a  finite  integration  time.  For  digital  signal  processing  it  is 
convenient  to  represent  the  process  by  its  discrete-frequency  samples.  This  leads  to  the  discrete 
Fourier  transform  (DFT).  The  DFT  pair  is  given  by 


x(«)  =  MN  X(k) 

*(")  «  ;  for  n,k  =  0,l,..,Ar-l.  /  ^ 

The  DFT  requires  complex  multiplication  operations  and  N(N-1)  complex  additions.  The 
fast  Fourier  transform  (FFT)  implements  the  DFT  with  fewer  multiplications.  The  FFT  has 
computation  complexity  (number  of  multiplications)  of  N/2  logj  N.  Time  or  frequency 
uncertainty  can  be  reduced  by  using  overtyping  techniques  [20]. 
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2.7  Time  Frequency  Distributions 

A  Time  Frequency  Distribution  (TFD)  [38]  can  describe  non-stationaiy  signals  by 
displaying  the  energy  density  as  a  function  of  time  and  frequency.  The  most  popular  TDF  is  the 
Wigner-Ville  Distribution  (WVD)  [15,38,39,40].  The  WVD  of  the  signal  s(t)  is  defined  by: 

TrVDg  (t,0))  =  I"  s*(t-T:/2)  s(t+t/2)  e dr;  (5) 

where  o)  =2  ic  f.  The  WVD  of  the  sum  of  two  signals  (i.e.,  s(t)  =  s,  (t)  +  Sj  (t))  is  given  by 
W(t,(o)  =  W,i  (t,Ci) )  +W22  (t,(o)  +W,2  (t,ti))  +W21  (t,(i>),  where 

it,w)  =  ^  J“s,\t-xl2)  Sj(t+xl2)  e  dx  (6) 

is  a  cross  term.  The  cross  term  is  complex-valued,  but  W,2  =  W2,*  and  hence 
W12  (1,0))  +  W2,  (t,(i))  is  real  and  W  (t,o))  =  Wj,  (t,o))  +W22  (t,o))  +2  Re  [W12  (t,o))  ],  where 
the  cross  term  2  Re  [W,2  (t,o))  ]  is  an  interference  term.  Window  functions  can  be  chosen  to 
improve  the  WVD  by  minimizing  cross  terms. 

2.8  Wavelets 

Wavelet  analysis  is  a  new  approach  to  represent  non-stationary  signals.  In  the  following 
discussion  we  introduce  wavelet  analysis  concepts.  The  function  f(t)  may  be  expressed  as 

M  =  E*  ^^(f)  ; 

with 

>  (7) 

=  f"  AO  v^(t)  dt  ; 

J  *>00 

where  file  set  Vk(t)  typically  forms  an  ortho-normal  set  and  <>  notation  denotes  an  inner 
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product,  i.e.,  a  projection  of  the  time  data  onto  the  k*  basis  function. 

There  are  different  t3/pes  of  wavelets;  orthogonal,  non-orthogonal,  and  bi-orthogonal 
wavelets.  The  Daubechies  family,  Symmlet,  Coiflet,  and  Meyer  wavelets  are  examples  of 
orthogonal  wavelets,  while  the  Morlet  wavelet  is  an  example  of  a  non-orthogonal  wavelet 
[16,34,35,37,42-47]. 

2.9  Continuous  Wavelet  Transform 

The  continuous  wavelet  transform  (CWT)  forms  the  mathematical  basis  for  wavelet 
analysis.  In  the  wavelet  analysis  all  basis  functions  can  be  generated  from  a  single  function 
called  the  mother  wavelet,  which  is  usually  denoted  by  i|t(t).  The  other  wavelets  can  be  generated 
using  two  distinct  operations;  scaling  and  translation.  Scaling  is  the  dilation  or  compression  of 
the  wavelet  function  according  to  a  specific  scaling  value.  The  scale  is  denoted  by  s.  The 
translation  allows  shifting  of  the  (scaled)  wavelet  to  a  desired  position  in  time.  This  shift  is 
denoted  by  a.  The  scaled  and  translated  wavelet  is  denoted  by 

where  Vs  is  a  nonnalization  factor.  The  integral  form  of  CWT  of  the  finite  energy  signal  f(t) 
with  respect  to  the  wavelet  function  t|r(t)  is  given  by  [42,46] 

W^(s,a)  = 

r'r-  (9) 

=  l/y/s  f  fit)  TK(t-a)/s)  du 

J  -oo 

The  wavelet  analysis  computes  inner  products  of  the  signal  and  the  wavelet  functions.  We 
can  also  interpret  the  wavelet  analysis  as  a  linear  operation  which  transforms  file  signal  nging 
modified  kernel  functions.  The  kernel  of  the  transform  is  file  mother  wavelet,  and  the 
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modifications  are  the  scaling  and  translation  operations  [44,45].  The  wavelet  operation 

can  be  interpreted  as  a  bandpass  fimction  which  implies  that  the  wavelet  must  be  an  oscillatory 
fimction. 

2.10  Scalogram 

Using  the  scalogram,  a  signal  f(t)  is  characterized  by  the  distribution  of  |  Wf  (s,a)p  over 
the  time-scale  plane.  The  quantity  |  Wj-  (s,a)p  may  be  viewed  as  a  spectral  density  in  units  of 
power  per  scale  [44].  Consequently,  the  scalogram  represents  the  power  spectral  density  of  the 
signal  over  the  time-scale  plane.  The  quantity 

r  da 

represents  the  portion  of  the  signal  energy  contained  within  the  scale  s.  Here 


L 


<fO)  <  oo 


and  T(g))  is  the  Fourier  transform  of  i|r(t).  This  fact  is  exploited  in  identifying  the  scale  for  each 
firequency  hop.  In  summary,  the  CWT  is  a  linear,  time-shift-mvariant,  time-scaling-invariant, 
and  firequency-scaling-invariant  operator. 
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2.11  Discrete  Wavelet  Transform 

The  CWT  is  defined  by  an  integral  transform  over  continuous  variables  in  the  scale  s  and 
the  time  shift  a.  In  practice  a  discrete  grid  for  s  and  a  is  used.  A  widely  accepted  discretization  is 
to  specify  s=sj^  and  k=na„s„'",  where  m  and  n  are  integers,  So>  1,  and  ao>  0  [45].  Furthermore, 
if  we  select  s„  =  2  and  a^  =  1  we  obtain  the  well-known  dyadic  wavelet  sampling  (tiling)  grid. 
Hence  the  scaled  and  translated  wavelet  indexed  by  m  and  n  is  given  by 

^m.„(t)  =  2-'^i|;(2-"t-n). 

We  will  briefly  introduce  the  WT  firom  the  perspective  of  the  multi-resolution  (MR)  analysis. 

The  signal  (or  the  time  function)  f(t)  is  expanded  in  terms  of  the  wavelet  functions.  These 
wavelets  have  a  firequency  bandpass  shape,  so  they  result  in  a  set  of  successive  details  of  the 
signal.  For  the  approximation  we  need  a  special  basis,  called  the  scaling  function  4)  (t),  which  is 
not  a  wavelet.  It  has  a  low  pass  firequency  behavior  and  performs  averaging.  The  discrete  wavelet 
synthesis  equation  is  given  by 

XO  =  Sr=-  «?(^)<i>i(0  +  5^"o  £r=o 

where]  and  k  are  integers,  the  coefficients  c(k)  constitute  the  coefficients  of  the  approximation, 
while  d(j,k)  constitutes  the  coefficients  of  the  added  details  or  equivalently  the  fine  resolutions 
[47].  If  the  wavelets  and  the  associated  scaling  functions  form  an  ortho-normal  set  of  basis 
functions  the  coefficients  are  given  by  c(k)  =  <  fi[t),4)k  (t)  >,  and  d(jdc)  =  <  f(t),4rjjj(t)  >. 

Here  4>(t)  is  a  lowpass  function  whose  firequency  response  is  the  same  as  the  fi'equency  response 
of  i|r(t)  excq)t  that  the  firequency  center  of  the  bandpass  filter  is  shifted  to  baseband  (i.e.,  centered 
at  DC). 
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2.12  Scaling  and  Wavelet  Equations 

The  Multi-resolution  subspace  representation  leads  to  a  method  to  formulate  two 
equations  in  terms  of  the  unknown  scaling  and  wavelet  functions: 

<l>(0  =  ^2  Ei  *o(*) 

This  equation  includes  two  different  scales  of  the  scaling  function  and  is  known  as  the  scaling  or 
dilation  equation.  The  coefficients  h„  (n)  are  called  the  scaling  filter  coefiBcients.  We  also  have 

m  =  E*  *i(*) 

which  is  called  the  wavelet  equation.  The  coefficients  hj  (n)  are  the  wavelet  filter  coefficients 

[42.47] . 

For  discrete  data  the  filter-bank  concept  leads  to  a  simple  method  for  computing  the 
wavelet  coefficients.  The  wavelet  function  is  replaced  by  the  coefficients  of  the  wavelet  filter 
hj(n)  and  the  scaling  function  is  replaced  by  the  coefficients  of  the  scaling  filter  ho(n).  In 

[42.46.47] ,  the  following  two  recursive  equations  are  obtained: 

=  E„  c(j+l,m) 

‘^(/»*)  =  Em  ^i(»*"2*)  c(/+1,ib)  . 

These  two  recursive  equations  enable  us  to  compute  the  j*  scale  wavelet  transform  and  is  known 
as  Mallat's  algorithm  [46]. 

2.13  Daubechies  Wavelet  Family 

The  Daubechies  wavelet  family  is  a  compactly  supported  ortho-normal  set  of  wavelet 
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functions  [42].  The  Daubechies  wavelet  are  obtained  by  solving  the  scaling  and  wavelet 
equations.  An  additional  set  of  constraints  is  applied  to  satisfy  the  maximum  number  of 
vanishing  moments  for  each  wavelet.  This  report,  and  the  Matlab  wavelet  toolbox  [59],  uses  the 
notation  that  a  Daubechies  wavelet  of  order  N  has  2N  coefficients.  The  wavelet  of  order  N  has 
finite  support  over  [0,2N-1],  or  equivalently  the  corresponding  FIR  filter  has  2N  multipliers.  The 
number  of  vanishing  moments  is  an  indication  of  the  smoothness  of  the  wavelet  filter.  The 
number  of  vanishing  moments  implies  the  number  of  zeros  of  T(o))  at  co  =  tt.  The  higher  the 
order  the  longer  and  smoother  the  Daubechies  filter  will  be. 
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3.  WAVELET  TRANSFORMS  AND  CORRELATION  FUNCTIONS 

Correlating  two  functions  provides  a  measure  of  their  similarity.  The  Wiener-Khinchin 
theorem  relates  the  signal's  auto-correlation  function  and  power  spectral  density  for  a  stationary 
process.  Wavelet  decomposition  can  be  used  to  represent  non-stationary  signals  over  the 
time-scale  plane.  We  will  examine  the  wavelet  transform  of  the  correlation  function  as  an 
alternative  for  non-stationary  signal  representation. 

3.1  Correlation  Functions 

Depending  on  the  underlying  process,  various  definitions  can  be  given  to  the 
auto-correlation  function  (ACF).  The  process  may  be  deterministic,  stochastic,  stationary  or 
non-stationary. 

The  ACF  of  a  stochastic  process  is  the  correlation  between  two  samples  of  the  process 
taken  at  t,  and  tj,  and  is  defined  as  R(t„t2)  =  E{x  (t,)  x*(t2)}  ,  where  E  {  }  is  the  expectation 
operator  and  *  stands  for  the  complex  conjugation.  For  a  stationary  (i.e.,  wide-sense  stationary) 
process,  R(t„t2)  depends  only  on  the  time  lag  t=  t,  -tj.  The  Wiener-Khinchin  theorem  defines 
the  relationship  between  the  correlation  function  and  spectral  density  as 
Sxx(“)=  J  R(t)  exp(-j  0)  t)  dr. 

3.2  The  Instantaneous  Correlation  Function 

The  ACF  of  a  deterministic  or  stochastic  process  is  computed  using  time  domain 
averaging  or  the  expectation  operator,  respectively.  This  means  that  a  smoothing  process  has  to 
be  applied  to  compute  the  correlation  functions.  The  instantaneous  correlation  function  (ICF) 
does  not  use  an  averaging  operation.  The  instantaneous  correlation  function  is  simply  defined  as 
the  product  of  two  samples  of  the  process.  These  two  samples  are  drawn  at  two  time  instants 
centered  about  time  t.  The  instantaneous  correlation  function  R(  t,T)  is  defined  as 
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R'(  t,T)  =  X  (t+T/2)  X*  ( t-  t/2),  where  i  stands  for  the  instantaneous  nature  of  the  correlation 
function  [49]. 

If  x(t)  is  a  sinusoidal  signal  then  the  multiplication  to  obtain  the  values  of  R*(  t,T) 
generates  cross  terms  in  the  ICF.  For  example,  the  real-valued  smusoidal  signal 
x(t)  =  A  cos  (o)  t)  has  an  ACF  given  by  R(t)  =  A^  /2  cos  (o)  t),  while  the  ICF  is  given  by 
R'(  t,T)  =  AV2  [cos  (2  (0  t’s)  +  cos  ( to  t)  ]. 

The  ACF  of  a  single  sinusoidal  signal  has  only  one  component  and  no  cross  term,  while  the  ICF 
has  cross  terms.  If  the  signal  x(t)  is  represented  by  its  analytic  form,  say  x(t)=A  exp  (j  to  t),  then 
its  ICF  is  given  by  R'(  t,T)  =  A^  exp  (j  o  t). 

That  is,  the  ICF  of  a  single  complex  exponential  signal  has  no  cross  term.  To  minimize  cross 
terms  from  the  negative  frequency  components  we  use  the  analytic  form  of  the  data  [49]. 

3.3  WAVELET  TRANSFORMS  OF  CORRELATION  FUNCTIONS 

The  wavelet  transform  of  the  stochastic  ACF  of  a  stationary  process  is  addressed  in  [48]. 
The  wavelet  transform  of  the  ACF  of  a  detaministic  signal  will  have  a  .similar  expression.  Let 
Wxx  (s,a)  denote  the  wavelet  transform  of  R(t).  Note,  the  subscript  in  W^x  (s,a)  stands  for  the 
wavelet  transform  of  the  ACF  of  x(t)  in  contrast  to  (s,a)  which  denotes  the  wavelet 
transform  of  x(t).  The  wavelet  basis  function  is  denoted  by  g(T).  The  wavelet  transformation 
will  transform  the  lag  variable  t  to  the  shift  variable  a  and  the  scale  s.  W^x  (s,a)  ,for  positive  s 
is  given  by: 

=  -rf” 

s 

=  ^  /.I 
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This  equation  has  the  form  of  an  inverse  Fourier  transform  from  the  variable  f  to  the  variable  a. 
We  can  write  W^x  (s,a)  =  {A  Sxx  (F)  G*(sf)  }  and  deduce  that  the  wavelet  transform,  at 
any  scale  s  >  0,  represents  a  linear  filtering  operation  using  a  band  pass  filter  whose  impulse 
response  is  the  (time-reversed)  wavelet  function  at  scale  s.  Equivalently,  the  filter  has  a 
frequency  response  given  by  the  FT  of  the  scaled  wavelet. 

The  wavelet  transform  of  the  ACF,  R(t  ),  of  the  stationary  finite-energy  signal  x(t), 
gives  a  band  pass  filtered  version  of  the  power  spectral  density  Sxx  (0  of  this  signal  (up  to  a 
constant,  v^s,  the  band  pass  filter  used  dependents  on  the  chosen  wavelet  function  and  scale. 

3.4  The  Wavelet  Transform  of  the  Instantaneous  Correlation  Function 

The  Wigner-Ville  Distribution  (WVD)  is  used  to  represent  non-stationary  processes. 

The  WVD  applies  a  one-dimensional  Fourier  transformation  to  the  ICF.  The  Fourier  transform 
takes  the  delay  x  to  the  frequency  f,  leaving  the  global  time  variable  t  unchanged.  This  allows 
the  display  of  the  time  evolution  of  the  spectrum  of  the  signal.  For  one-dimensional  time  signals, 
the  one-dimensional  wavelet  transform  carries  out  a  transformation  from  one  global  time  variable 
t  to  the  two  wavelet  variables,  the  shift  a  and  the  scale  s.  Consequently,  the  signal  is  represented 
by  a  time-scale  distribution  in  the  wavelet  domain.  The  wavelet  domain  is  called  the  time-scale 
domain.  For  the  two-dimensional  surface,  indexed  by  time  t  and  the  delay  x,  we  carry  out  the 
wavelet  transformation  along  the  delay  axis.  This  permits  a  display  as  a  function  of  time.  Let 
Vx  (t,f)  denote  the  WVD  of  the  signal  x(t), 

'^x  =  f  x(it-xl2)  X  *(r+x/2)  e  dx 

=  dx  , 
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and  let  the  Fourier  transform  of  the  wavelet  basis  function  be  given  by 

=  s  G(5f) 
s 

so 

gi—)  =  r  s  G(sJ)  df, 

s  j  -m 

then  using 

^  XX  ~  ~~Z  f  g  (  )  » 

we  have 

W  f  G  Us/)  (tj)  df . 

This  equation  is  in  the  form  of  an  inverse  Fourier  transform.  ' (t;  s,a)  and 

fs  G*  (sf)  (t,0  are  a  Fourier  transform  pair  with  respect  to  the  variable  a  and  f.  This  relation 
suggests  that  we  can  obtain  a  filtered  version  of  fiie  WVD  by  Fourier  transforming  '  (t;  s,a). 

3.5  FREQUENCY  HOPPED  SIGNALS  AND  THEIR  CORRELATION  FUNCTIONS 
Communication  systems  can  utilize  a  large  number  of  digital  modulation  techniques; 
spread  spectrum  modulation  being  one  of  them.  Spread  spectrum  refers  to  any  modulation 
scheme  that  produces  a  transmitted  bandwidth  much  larger  than  the  information  bandwidth  [50]. 
We  will  briefly  address  different  digital  modulation  schemes  and  focus  on  firequency  hopping. 
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3.6  SPREAD  SPECTRUM  COMMUNICATION  SIGNALS 

Spread  spectrum  (SS)  commumcation  signals  are  characterized  by  a  wide  transmission 
bandwidth  and  a  low  power  spectral  density  [19,  50,51].  SS  signals  have  two  main  advantages: 

i)  The  message  has  a  low  probability  of  being  intercepted  (LPI)  as  a  result  of  the  wide  frequency 
band  and  the  low  power  spectral  density  of  the  signals. 

ii)  .SS  systems  can  reject  jamming  signals  and  allow  users  to  share  the  same  frequency  band. 
Among  the  different  possible  SS  modulation  formats,  the  following  three  are  prevalent: 
Frequency  Hopping  (FH),  Direct  Sequence  (DS)  Modulation  and  Time  Hopping  (TH). 

3.7  THE  INSTANTANEOUS  CORRELATION  FUNCTION  OF  FREQUENCY 
HOPPED  SIGNALS 

Spread  spectrum  studies  usually  consider  the  FH  signal  as  a  stationary  process 
[19, 50,64]  even  though  the  spectrum  of  the  FH  signal  varies  with  each  hop  interval.  The 
stationary  correlation  representation,  using  time  averaging,  is  not  suitable  for  this  process. 

One  way  to  identify  the  FH  signal  is  to  monitor  the  time-frequency  evolution  of  the 
signal.  Hence,  we  need  to  keep  the  time  dependency  in  flie  correlation  representation.  This  is 
achieved  by  using  the  instantaneous  correlation  function  (ICF). 

The  FH  signal  may  be  represented  as  successive  intervals  (i.e.,  hops)  of  single  frequency 
complex  erqronential  (Fig.  2a).  The  frequency  within  each  interval  (i.e.,  flie  hop  frequency)  is 
controlled  by  a  random  (but  known  to  the  user)  sequence.  We  assume  without  loss  of  generality  • 
that  any  two  successive  hops  will  have  different  frequencies. 

The  frequency  difference  of  adjacent  hops  will  generate  the  patterns  of  the  instantaneous 
correlation  functions.  The  IFC  for  values  of  |t1  ^  T^  (i.e.,  allows  correlating  adjacent  hops 
only)  and  using  values  of  (t-i-T/2)  and  (t-T/2)  such  that  the  values  are  confined  to  be  within  the 
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Figure  2a:  Time  behavior  of  an  FH  signal. 
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Figure  2b:  FH  signal  and  the  cellular  (diamond)  stmcture  of  the  ICF. 
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L*  hop,  is  given  by;  R(t,T)  =  e^‘”, 

where  o)  is  the  radian  frequency  of  the  complex  valued  sinusoid.  We  note  that  the  values  of 
(t+T/2)  and  (t-T/2)  are  confined  to  be  within  the  same  L*  hop  if  they  satisfy 
(L-l)Th  <  (t+T/2)  and  (t-T/2)  <  LTj, .  This  inequality  forms  the  boundaries  of  a  diamond 
pattern  for  a  given  value  of  L.  This  cellular  structure  is  shown  in  Figure  2b.  Inside  each  diamond 
the  ICF  is  obtained  by  correlating  signals  from  the  same  hop,  while  outside  the  diamond  the  ICF 
is  obtained  by  correlating  signals  belonging  to  two  consecutive  hops. 

The  correlation  function  R„  „  (t,T)  (i.e.,  in  region  )  is  given  by: 

Rm,n  =  exp  j  ((i)„  -(oj  t  +(ci)„  +  (dJ  t/2;  where  m  and  n  are  the  indices  of  the  two  adjacent 
hops.  We  note  that  within  the  main  diamond  of  the  n*  hop,  (i.e.,  m=n),  the  ICF  is  given  by 
Rm,n  =  exp  j  ((•)„  t)  while  outside  the  main  diamond  (i.e.,  in  the  upper  triangle  between  hop 
numbers  n  and  n+1),  the  ICF  is  given  by  (t,T)  =  exp  (j  (w^,  -wj  t  +((o^,  +  t/2}  . 

The  lower  half  of  the  ICF  has  hermitian  symmetry  relative  to  the  upper  half  and  does  not  provide 
any  additional  information. 

In  summary,  the  instantaneous  correlation  function  of  a  FH  signal  will  exhibit  cellular 
(i.e.,  diamond)  patterns.  Inside  the  L*^  diamond  the  ICF  has  a  single  complex  Exponential 
component  along  the  delay  axis  representing  the  L*  hop  frequency.  Outside  flie  diamond,  RXt,T) 
is  a  product  of  two  terms,  exp  j  (co„  -o)  J  t  and  exp  j  (o)„  +  w  J  t/2  ,  where  o)„  and  <o„  are  two 
consecutive  hop  frequencies. 
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4.  PROCESSING  SCHEME 

The  wavelet  transform  generates  one  surface  for  each  scale.  One  can  visually  inspect  the 
wavelet  surfaces  to  identify  the  FH  signal  and  obtain  an  estimate  for  the  hop  time  interval. 
Alternatively,  one  can  use  a  processing  scheme  to  estimate  the  hop  start/stop  times,  the  hop-scale 
pattern,  and  the  hop  frequency.  For  the  extraction  of  the  hop  start/stop  times  an  edge  detector  is 
used.  An  estimate  of  the  hop-scale  pattern  can  be  obtained  by  performing  an  energy  analysis. 

The  energy  analysis  assigns  a  scale  index  (called  the  proper  scale)  to  each  hop.  The  proper 
scale,  is  that  scale  which  has  the  greatest  energy  content  (i.e.,  spectral  components  live  in  the 
spectral  region  covered  by  the  scale  under  consideration).  The  sequence  of  proper  scales, 
representing  the  hop  sequence,  is  called  the  hop-scale  pattern. 

If  a  hop-scale  pattOTi  is  detected,  it  provides  the  evidence  that  a  frequency  hopped  signal 
is  present.  If  some  or  all  frequencies  of  an  FH  signal  reside  in  the  spectral  region  of  one  scale 
then  a  follow-on  spectral  estimation  will  indicate  different  frequency  components  as  a  function 
of  the  hop  intervals  which  still  permits  the  identification  of  the  FH  process. 

For  FH  signals,  the  ICF  displays  a  cellular  pattern  (i.e.,  see  Fig.2),  where  each  hop 
results  in  a  diamond  pattern  with  a  width  equal  to  the  hop  interval.  The  diamond  intersects  the 
time  axis  at  the  hop  start/stop  times.  The  wavelet  surface  at  the  proper  scale  emphasizes  signals 
which  belong  to  that  scale,  while  other  signals  (i.e.,  out  of  band  components)  are  attenuated. 

The  interception  problem  usually  assumes  some  prior  knowledge  about  the  signal  of 
interest.  In  our  case,  we  assume  the  minimum  and  maximiun  hopping  frequencies  are 
approximately  known  and  the  data  is  propwly  sampled. 
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4.1  Discrete-Time  Implementation  of  the  Instantaneous  Correlation  Function 

Let  R(n,u)  define  the  ICF  of  the  discrete-time  signal  x(n),  given  by: 

R(n,u)  —  X  (n+u/2)  x*  (  n-u/2),  where  n  is  time  and  u  is  the  delay.  The  combined  index 
n±  u/2  should  assume  integer  values.  We  can  insert  zeros  into  the  ICF(i.e.,  R(n,u),  for  odd  u, 
is  set  to  zero),  or  we  can  let  u  =  2m.  The  last  approach  is  the  one  adopted  in  this  report. 

A  one-dimensional  wavelet  transform  is  performed  in  the  direction  of  the  delay  u  for 
each  time  element  of  the  correlation  function.  The  Matlab  wavelet  toolbox  [59]  is  used  with  the 
convention  that  the  highest  band  of  passband  frequencies  is  denoted  by  scale  number  1. 

4.2  VISUAL  IDENTIFICATION 

We  investigated  different  types  of  wavelets  as  well  as  different  surface  representations.  A 
complex-valued  Daubechies  wavelet  of  order  3  is  given  in  [60]  and  is  used  for  comparison  with 
the  real  valued  Daubechies  wavelet  of  the  same  order.  Operating  on  the  complex  valued  ICF  with 
a  real  or  complex  valued  wavelet  results  into  a  complex  valued  scale  (surface)  output.  A  complex 
valued  surface  can  be  represented  by  its  magnitude,  phase  or  potentially  the  real  or  imaginary 
component.  A  visual  identification  technique  was  used  to  select  the  type  of  wavelet  and  the  type 
of  scale  representation.  Based  on  a  large  number  of  simulations  in  conjunction  with  an  opinion 
test  it  was  concluded  diat  a  real  valued  wavelet  will  suffice  and  the  surfaces,  using  the  real 
output,  provide  superior  results.  Hence  in  all  follow-on  work,  real  valued  wavelets  and  real 
valued  scale  surfaces  are  used.  This  processing  approach  was  also  verified  using  the  Morlet 
wavelet  [56]. 

The  objective  of  the  opinion  test  is  to  identify  the  frequency  hop,  interpreting  the  cellular 
structure  in  the  wavelet  surface,  and  to  identify  the  diamond’s  time  axis  intercepts  to  serve  as  an 
estimate  of  the  hop  start/stop  times. 
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Figures  4.1  and  4.2  show  examples  of  the  real  part  of  the  WT  obtained  with  the 
real-valued  Daubechies  wavelet  of  order  3.  Figure  4.1  and  4.2  display  the  real  part  of  the  IFC  and 
the  real  part  of  the  first  5  scales  of  a  FH  signal  with  an  SNR  of  10  and  3  dB,  respectively.  The 
frequencies  of  the  FH  signal  are  such  that  they  hop  the  scales  in  a  staircase  fashion  (i.e.  time 
segment  1  or  hop  1  is  on  scale  1,  time  segment  2  or  hop  2  is  on  scale  2,  etc,).  We  note  that  we  are 
unable  to  identify  the  FH  structure  from  the  original  ICF  surface,  denoted  by  "CF.”  Wavelet 
surfaces,  labeled  "Sir", ..,  "S5r",  allow  identification  of  diamond  patterns  at  hop  number  1, 5 
,  respectively.  The  diamond  patterns  are  detectable  since  they  are  presented  by  contour  lines  of 
constant  heigjit  that  run  from  left  to  right.  If  we  were  to  plot  the  delay  trace  that  goes  through  a 
diamond,  for  a  given  fixed  value  of  time,  we  will  observe  a  sinusoidal  pattern  having  a  fixed 
fi-equency.  These  figures  demonstrate  that  we  can  identify  the  FH  structure  fi:om  the  wavelet 
surfaces,  while  it  is  not  possible  to  do  so  firom  the  ICF  surface.  At  high  SNR’s  (i.e.,  10  dB  or 
better),  we  can  also  determine  the  hop  tunes  (i.e.,  where  the  diamond  intercepts  the  time  axis) 
easily.  A  color  coded  display  is  superior  to  the  black  and  white  representation  used  in  this  report. 
4.3  ENERGY  ANALYSIS  AND  SCALE  IDENTIFICATION 

If  correct  hop  timing  information  is  available,  we  can  perform  an  energy  analysis  for  each 
hop.  Parseval’s  theorem  for  the  complete  orthogonal  filter  bank  (over  L  partitions)  is  applicable 
to  the  discrete  time  wavelet  analysis  [45].  So  ||x(n)|p=  (  |C(L,2a)p+  Sj.,  ^  |d0’,2a+l)p  ), 

where  C(L,2a)  are  the  scaling  coefficients  at  the  scale  L,  d(j,2a+l)  are  the  wavelet  coefficients  at 
scale  j,  and  a  is  the  wavelet  shift  variable.  The  quantity  E(j)  =  S  "  |  d(j,a)  |  ^  represents  the 
signal  energy  over  the  j*  scale.  Energy  per  sample  is  defined  as:  A(j)  =  EQ)/  N(j),  whore  j  is  the 
scale  index,  E(j)  is  the  total  energy  at  the  j*  scale,  and  N(j)  is  the  number  of  wavelet  coefficients 
at  this  scale. 
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Figure  4.1  Wavelet  Surfaces  of  FH  Signals  Using  Daub-3  at  10  dB  (scale  1-5). 
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To  identify  the  proper  scale  one  can  use  the  maximum  value  of  the  wavelet  coefficients, 
the  total  energy,  or  the  energy  per  sample.  Table  1  summarizes  the  energy  distribution  obtained 
via  wavelet  analysis  using  Daubechies  wavelet  of  order  2  (Daub-2)  and  order  10  (Daub- 10).  The 
signal  is  64  samples  long  and  consists  of  three  sinusoids.  The  first  signal  has  a  frequency  of  3/8 
Fj,  the  second  has  3/16  Fj,  and  the  third  has  3/32  Fj,  where  Fj  is  the  sampling  frequency.  Hence 
the  first,  second,  and  third  input  signal  is  contained  in  the  first,  second,  and  third  scale, 
respectively.  Using  the  information  from  table  1  we  can  conclude  that  the  total  energy  of  the 
input  signals  is  distributed  across  the  scales  and  that  the  sum  of  the  total  energies  over  the  scales 
is  slightly  less  than  the  total  signal  energy  since  we  disregard  any  contribution  from  the  low  pass 
section.  The  proper  scale  (where  the  signal  resides)  has  the  greatest  share  of  the  total  energy. 

This  share  increases  with  an  increase  of  the  wavelet  length.  This  is  attributable  to  the  smaller 
spectral  leakage  of  longer  wavelets.  Energy  per  sample  at  the  proper  scale  is  larger  if  the  signal 
resides  at  a  higher  scale.  The  gain  factor  in  energy  per  sample  (if  the  signal  resides  within  scale  2 
rather  than  scale  1)  is  about  1.41  and  1.64  for  Daub-2  and  Daub-10,  respectively.  Ideally,  the 
gain  factor  in  energy  per  sample  should  be  2  per  scale  index.  For  an  automated  detection  scheme 
the  energy  or  energy  per  sample  is  the  superior  test  statistics.  To  declare  that  a  signal  belongs  to  a 
particular  scale,  that  scale  must  be  the  dominant  one  relative  to  all  otihier  scales;  The  larger  the 
ratio  the  better  the  discrimination.  We  note  in  table  1  that  the  energy  and  energy  per  sample 
outperform  the  maximum  coefficient  at  all  scales  and  all  wavelet  sizes  tested. 

431  Hop-Scale  Pattern 

The  scale  identification  assign?  a  scale  index  to  each  hop  of  the  observed  signal.  The  hop 
is  assigned  to  the  scale  whose  energy  per  sample  is  the  largest  among  the  values  of  the  other 
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Table  1:  Energy  Distribution  of  Daub-2  and  Daub-10  Wavelets 


Measure 

Wavelet 

Scale 

II 

00 

Sign^  with 
|/  =  (3/16)F, 

/  =  (3/32)F, 

1 

1.3365 

Daub-2 

.2 

0.5753 

Max 

3 

0.3743 

1.8480 

Coefficient 

1 

1.2524 

Daub-10 

2 

3 

0.1854 

2.8296 

r?- 

1 

29.7687 

7.2565 

Daub-2 

2 

1.1313 

23.3282 

8.2775 

Total 

3 

1.0495 

0.9492 

21.2601 

Energy 

1 

31.4508 

1.6228 

Daub-10 

2 

0.3980 

29.9017 

1.9671 

3 

0.0875 

0.3900 

29.0447 

1 

Sample 

Daub-2 

2 

Energy 

3 

0.1049 

2.1260 

(average 

1 

0.7671 

per  sample) 

Daub-10 

2 

0.0133 

3 
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scales.  This  should  correspond  to  the  true  scale  location  of  the  frequency  hop  and  is  called  the 
proper  scale.  A  sequence  of  hops  will  result  in  a  sequence  of  proper  scales  forming  the  hop-scale 
pattern. 

4.32  Success  Rate 

The  performance  of  scale  identifrcation  is  evaluated  via  the  success  rate  Pj^.  This  is  done 
by  generating  known  hop-scale  patterns  and  obtain  the  percentage  of  foe  correctly  identified 
hop-scales.  The  success  rate  is  defined  as;  Py  =  number  of  correct  hop-scales  /  total  number  of 
hops.  The  quality  of  scale  identification  depends  on  foe  height  of  foe  greatest  energy  per  sample 
relative  to  foe  other  sample  energies  from,  other  scales. 

The  spectral  density  of  foe  ICF  for  a  white  noise  input,  in  foe  delay  direction,  has  a 
trian^ar  spectral  shape.  To  allow  comparison  foe  spectral  shape  foe  energy  per  sample  for  all 
wavelet  scales  must  be  corrected  accordingly.  Performance  in  terms  of  foe  success  rate  Py  is 
given  in  chapter  5. 

4.4  FREQUENCY  ESTIMATION 

The  Fourier  Transform  of  foe  wavelet  surface  gives  a  bandpass  filtered  version  of  foe 
WVD.  The  hop  frequency  can  be  obtained  from  a  Fourier  transform  over  foe  delay  direction  of 
foe  wavelet  surface  (i.e.,  over  foe  main  diagonal  of  foe  appropriate  diamond  region). 

The  Fomier  transform  of  foe  wavelet  coefficients  can  be  used  as  a  spectral  estimate  (i.e., 
petiodogram)  at  all  scales.  The  frequency  resolution  dependents  on  foe  parameters  of  foe  Fourier 
transform. 

The  frequency  resolution  (i.e.,  foe  minimum  spacing  between  two  resolved  narrow  band 
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components)  of  the  DFT  is  approximately  equal  to  Af  =/N.  At  any  given  scale  k,  the  number  of 
data  points  N(k)  is  related  to  the  number  of  input  data  points  N  by  N(k)  =  N/  2'‘.  The  sampling 
frequency  of  the  scale  output  (i.e.,  detail  function)  is  scale  dependent,  i.e.,  (k)  =  ,  where 

Fj  is  the  input  sampling  frequency.  At  the  k*  scale,  both  the  number  of  data  points  (wavelet 
coefficients)  and  the  sampling  frequency  have  been  reduced  by  the  same  factor.  Consequently, 
the  FT  of  the  N(k)  data  points  has  a  frequency  resolution  that  is  constant,  independent  of  the 
wavelet  scale  being  addressed. 

4.41  Success  Rate 

Performance  of  the  frequency  estimation  procedure  is  evaluated  in  terms  of  the  success  rate,  Pf. 
The  hop  frequency  is  considered  correctly  estimated  if  the  spectral  peak  is  at  the  true  spectral  bin. 
Pf  is  defined  as:  Pf  =  number  of  correct  hop  frequencies  /  total  number  of  hops. 

The  hop  frequency  is  given  by  the  bin  number  corresponding  to  the  peak  of  the  magnitude  of  the 
FT  over  a  specified  region  of  the  wavelet  surface.  The  quality  of  the  frequency  estimation 
depends  on  the  spectral  height  of  the  peak  relative  to  the  average  backgrmmd. 

Given  the  correct  hop  start/stop  times,  the  hop  frequency  may  also  be  estimated  directly 
from  the  time  signal  or  from  the  ICF  surface.  To  extract  the  frequency  from  the  original  time 
signal  we  can  use  the  FFT  of  the  time  data  over  the  hop  length.  The  FFT  is  a  matched  filter  for 
sinusoidal  signals  in  white  Gaussian  noise.  Thus,  an  optimal  performance  is  expected  relative  to 
the  nonlinear  processing  of  the  signal  throu^  the  ICF  computation  and  the  linear  wavelet 
transformation. 

WensePf^and  P^  to  denote  the  success  rates  of  the  frequency  estimator  when  using  the 
original  time  signal  and  the  ICF  surface,  respectively.  Results  are  provided  in  chapter  5. 


28 


4.5  ESTIMATION  OF  HOP  TIMES 

We  recall  that  the  ICF  surface  and  wavelet  surfaces  have  a  cellular  structure  consisting  of 
diamonds,  where  each  diamond  is  associated  with  a  specific  hop.  The  diamond's  interception  of 
the  time  axis  defines  the  start/stop  point  of  the  hops.  The  diamond  width  corresponds  to  the  hop 
interval  T^.  The  sides  (edges)  of  the  diamonds  of  hops  are  mutually  parallel  and  spaced  by  the 
hop  interval  T^.  There  are  many  approaches  one  can  use  to  solve  the  problem  of  hop  time 
estimation.  In  what  follows  we  will  use  a  technique  based  on  an  edge  detection  operator. 

Edge  detection  is  a  fimdamental  problem  in  image  analysis  since  edges  help  in 
identifying  objects.  There  are  two  basic  types  of  edge  operators,  the  gradient  operators  and  the 
compass  gradient  operators  [62,63].  The  gradient  operator  measures  the  gradient  of  the 
two-dimensional  image  in  two  orthogonal  directions.  It  is  usually  applied  to  detect  edges  with 
unknown  directions.  The  compass  operator  measures  the  gradient  of  the  two-dimensional  image 
in  a  specific  direction  (i.e.,  ±  ti/4  in  the  ICF  or  scale  outputs). 

The  compass  operator  is  applied  to  the  upper  half  of  the  wavelet  surfaces  summing  up  all 
contributions  according  to  the  compass  weights.  The  maximum  value  is  extracted  and 
determines  the  point  where  the  compass  array  matches  an  edge.  To  make  the  data  applicable  to 
compass  operations  one  needs  to  add  a  positive  number,  equal  in  magnitude,  to  the  smallest 
(i.e.,  the  most  negative)  surface  value.  Results  are  presented  in  chapter  5. 
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5.  SIMULATIONS  AND  RESULTS 


This  chapter  provides  simulation  results  of  the  techniques  introduced  in  chapter  4. 
Chapter  5  deals  with  wavelet  transforms  exclusively  applied  along  the  delay  axis.  We  will 
address  visual  inspection,  scale  identification,  frequency  and  hop  time  estimation.  Results  for 
wavelet  processing  along  the  delay  axis  are  given  in  chapter  6. 

5.1  VISUAL  INSPECTION 

To  detect  and  to  identify  the  FH  modulation  we  performed  an  opinion  test  by  examining 
the  wavelet  surfaces  visually.  Ten  participants  were  involved,  each  one  was  asked  to  identify  the 
diamond  patterns  of  the  FH  signal  from  the  wavelet  surfaces  at  all  pertinent  scales  and  for  all 
hops.  Two  types  of  wavelets  were  used;  the  Morlet  wavelet  and  die  Daubechies  wavelet  of  order 
3.  Both  wavelet  types  were  used  in  their  real  and  complex  form.  Four  SNR  values  were  used; 

10, 6, 3  and  0  dB.  Four  different  surface  representations  were  examined;  the  real  part,  imagiaary 
part,  magnitude,  and  phase.  To  minimize  biasing  of  the  test  results,  all  participants  started  to 
identify  the  surfaces  going  from  the  lowest  SNR  to  the  highest  SNR  value.  More  detailed 
scoring  tables  and  the  scoring  code  are  given  in  [56  ].  The  FH  signal  occupies  the  first  five  scales 
at  different  hop  times. 

Table  2  shows  scoring  results  based  on  the  real  values  of  the  scale  surfaces  for  an  SNR  of 
10  and  3  dB.  10  dB  is  the  highest  SNR  value  tested,  while  3  dB  is  the  minimum  value  that  still 
provides  an  accqjtable  identification  score.  The  values  of  the  ratings  range  from  0.2  to  1.  Here  1 
indicates  perfect  identification  of  die  hop  diamond  patterns  at  their  proper  time  locations  while 
0.2  indicates  just  a  detection  of  a  hop  pattern  in  the  backgroimd  noise. 

The  visual  opinion  test  indicates  that: 
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Table  2:  Summaiy  of  the  Identification  Score  Using  the  Real  Part  of  the  Scales. 


Wavelet 

SNR 

Scales 

Type 

[dB] 

SI 

S2 

S3 

S4 

S5 

Real 

1.0 

1.0 

0.95 

M 

Morlet 

3 

1.0 

1.0 

0.95 

B 

Complex 

1.0 

1.0 

^9 

0.75 

Morlet 

3 

0.7 

0.8 

0.4 

Real 

B 

n 

1.0 

0.8 

B 

0.2 

Daub-3 

B 

1 

0.7 

i 

0.2 

Complex 

B 

1.0 

1.0 

1.0 

1 

Daub-3 

B 

0.7 

i 

H 

Scores  vary  between  0.2  and  1,  where 

1  :  perfect  identification  of  hop  diamonds. 

0.2:  just  distinction  of  hops  from  background  noise. 


i)  The  FH  signal  can  be  identified  from  the  wavelet  surfaces  by  its  cellular  structure  which  is 
dominant  at  the  proper  scales.  That  is,  each  scale  will  emphasize  the  hops  which  belong  to  the 
scale  and  attenuate  other  out  of  band  spectral  components.  The  (diamond)  cellular  structure  can 
be  used  in  the  visual  estimation  of  the  hop  start/stop  times. 

ii)  The  FH  signal  can  be  easily  identified  at  SNR  levels  of -3  dB  and  above. 

iii)  The  real  part  (or  imaginary  part)  provides  the  best  representation  for  visual  inspection. 

iv)  The  real  value  of  the  wavelet  function  provides  a  better  surface  representation  than  the 
complex  valued  wavelets.  We  tested  the  Morlet  and  Daubechies  (Daub-3)  wavelets  [56].  Other 
types  of  wavelets  may  perform  differently  but  were  not  evaluated. 

v)  Other  modulation  schemes  such  as  ASK,  PSK,  MFSK  and  noise  only  patterns  will  have 
patterns  residing  at  one  scale  only  or  have  no  discernible  patterns  at  all.  Typical  plots  are  given  in 
[56]  verifying  that  if  diamond  patterns  are  noted  at  different  scales,  FH  signals  are  present. 

5.2  SCALE  IDENTIFICATION 

A  processing  scheme  is  used  to  extract  hop  start/stop  times,  the  hop-scale  pattern,  and 
the  hop  frequency.  Initially,  we  investigate  scale  identification  and  hop  frequency  estimation 
assuming  that  correct  hop  timing  is  available.  Estimation  of  hop  start/stop  times  is  examined  in 
section  5.4.  For  section  5.2  and  5.3  the  parameters  of  the  simulations  are  as  follows: 

i)  Signal  pattern  length:  5  hops. 

ii)  Wavelet  scales  occupied:  SI,  S2,  S3,  and  S4  (i.e.,  one  scale  is  used  twice). 

iii)  Wavelet  types:  Daub-2,  Daub-4,  and  Daub-8. 

iv)  SNR  range:  -10  to  10  dB. 

v)  Number  of  realizations:  10  per  scale  per  SNR  value. 
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The  hop  frequencies  of  the  FH  signal  are  spaced  to  generate  the  hops  according  to  a 
known  fixed  scale  test  pattern.  The  test  pattern  occupies  scales  SI  through  S4.  The  wavelet 
surfaces  are  generated  from  the  ICF  surfaces  at  the  relevant  scales  (i.e.,  scales  1  through  4).  The 
total  energy  of  each  hop  at  each  scale  is  computed  and  the  energy  per  samples  is  obtained  by 
dividing  the  total  energy  by  the  number  of  wavelet  coefficients  at  each  scale. 

For  each  hop,  the  scale  with  the  greatest  energy  per  sample  is  designated  as  the  proper 
scale.  The  resultant  hop  pattern  is  compared  to  the  known  hop  pattern  and  the  probability  of 
correct  identification  is  computed.  To  avoid  bias  from  the  colored  noise  of  the  ICF  surfaces  an 
equalization  is  performed  at  all  scales  prior  to  the  estimation  procedure. 

Figures  5.1  to  5.3  show  the  performance  of  scale  identification  using  the  success  rate  Pj^. 
Results  are  obtained  for  Daubechies  wavelets  of  order  2, 4,  and  8  for  scales  SI  through  S4.  For 
the  scale  identification  performance  we  consider  the  minimum  SNR  at  which  is  still  unity  as 
the  figure  of  merit.  Over  all  tested  scales  the  success  rate,  Py  assumes  the  value  of  1  at  different 
minimum  input  SNR  values.  This  is  a  fimction  of  the  order  of  the  wavelet  and  the  scale.  Figure 
5.1  shows  that  the  performance  of  Pjj  obtained  from  Daub-2  achieves  a  Pj^  of  1  at  an  SNR  of 
-1  dB  at  all  scales,  hence,  -1  dB  is  considered  the  minimum  SNR  value  for  Daub-2.  The 
minimum  SNR  value  for  Daub-4  is  -2  dB  at  all  scales  as  shown  in  Figures  5.2.  Figure  5.3 
indicates  that  an  SNR  level  of -IdB  or  better  is  required  to  guarantee  a  P^  level  of  1 .  For  a  Py  of 
0.9  we  need  -1  dB,  -  2  dB,  and  -  3  dB  for  Daub-2,  Daub-4,  and  Daub-8,  respectively.  This  shows 
that  longer  wavelets  perform  better  than  shorter  ones  in  terms  of  Pid.  The  exception  of  that  case 
is  the  performance  at  scale  SI  to  S4.  The  performance  degradation,  as  the  length  of  the  wavelet  is 
increased,  may  be  due  to  a  non-ideal  equalization  of  the  ICF  ^ectral  shape. 
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1 


Figure  5.1  for  Daub-2. 
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Figure  5.2  Py  for  Daub-4. 
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Theoretically,  wavelet  surfaces  at  higher  scales  have  higher  SNR  values  than  those  at 
lower  scales  due  to  the  reduced  passband  regions  at  higher  scales.  Generally,  the  scale 
identification  performance  increases  with  increases  in  scale  number,  wavelet  length,  and  SNR. 
There  are  some  small  inconsistencies,  but  we  attribute  the  anomalies  to  the  non-ideal 
equalization  and  the  non-ideal  filter  transfer  fimctions  (in  a  spectral  sense). 

5.3  FREQUENCY  ESTIMATION 

We  earned  out  simulations  to  evaluate  the  performance  of  frequency  estimation  using 
the  data  specified  in  section  5.2.  The  hop  frequencies  are  estimated  by  taVing  the  FT  of  the 
wavelet  coefficients  located  at  the  center  of  the  diamond  patterns  in  the  direction  of  the  delay. 
The  bin  corresponding  to  the  peak  value  represents  the  estimated  hop  frequency.  The  estimated 
hop  frequency  is  compared  to  the  trae  hop  frequency  and  the  probability  of  correct  frequency 
extraction  (the  success  rate  Pf )  is  computed.  The  estimated  frequency  is  considered  correct  if  the 
estimation  error  is  less,  in  percentage  of  the  true  frequency,  than  1/N,  where  N  is  the  length  of 
the  vector  of  the  wavelet  coefficients.  Figures  5.4  to  5.6  plot  the  success  rate  Pf  obtained  for 
different  wavelets  as  a  function  of  scale.  For  the  frequency  estimation  we  consider  the  minimiiTn 
SNR  for  a  given  Pf  value  as  the  figure  of  merit. 

Figure  5.4  shows  that  the  Pf  value  for  Daub-2  is  1  at  an  SNR  equal  to  0  dB  at 
most  of  the  scales.  A  Pf  of  0.9  is  obtained  for  all  scales  for  SNR  levels  greater  or  equal  to  -1  dB. 
The  minimum  SNR  value  for  Daub-4  and  Daub-8  at  a  Pf  of  0.9  is  -2  dB,  as  shown  in  Figures  5.5 
and  5.6.  For  Daub-2,  Daub-4  and  Daub-8  at  a  Pf  of  1 .0,  the  miniTmun  SNR  level  is  0  dB  at  scales 
SI,  S2,  and  S3. 
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Figure  5.6  Pf  for  Daub-8. 
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Values  of  the  success  rate  Pf  show  that  the  hop  frequency  can  be  reliably  estimated  from  most 
wavelet  surfaces  at  an  SNR  s  0  dB  by  using  only  one  FFT  at  the  center  of  the  diamond  area  in 
the  direction  of  the  delay. 

Hopping  frequencies  may  also  be  estimated  directly  from  the  original  signal  or  from  the 
ICF.  Figure  5.7  plots  the  performance  and  Pf,  obtained  from  the  time  signal  and  from  the 
IGF,  respectively.  The  plots  are  indexed  by  the  SNR  and  assume  that  exact  estimates  of  hopping 
start/stop  times  are  available.  By  contrast,  using  the  wavelet  surface,  under  the  best 
circumstances  (i.e.,  scale  3,  Daub-8),  we  need  an  SNR  ^  -3  dB  to  obtain  perfect  performance. 
The  SNR  should  be  about  -3  to  -5  dB  for  using  the  raw  time  signal  (i.e.,  need  hop  timing 
information).  The  frequency  estimation  success  rate  using  the  ICF,  at  a  P^  of  1,  requires  an  SNR 
value  of  0  dB  or  better. 

This  shows  that,  assuming  exact  estimates  of  hopping  start/stop  times  are  available,  hop 
frequencies  may  be  estimated  by  processing  the  original  signal  at  lower  SNR  values  than  can  be 
achieved  using  the  wavelet  or  the  ICF  surfaces.  The  benefit  obtained  by  analyzing  the  ICF 
surface  by  wavelet  analysis  is  significant  in  case  of  unknown  hop  start/stop  times.  We  also  note 
that  fewCT  computations  are  needed  when  we  estimate  the  frequency  from  wavelet  surfaces.  This 
is  due  to  applying  only  one  FT  per  hop,  using  few  coefficients  (i.e.,  decimated  wavelet  ou^ut). 
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5.4  HOP  TIMES  ESTIMATION 


This  section  presents  hopping  time  estimation  results  obtained  using  the  compass 
operator  referred  to  earher.  The  wavelet  surface  can  be  represented  by  its  upper  half  plane.  Then 
the  areas  of  interest  (i.e..  Fig.  2)  have  a  triangular  pattern  instead  of  a  diamond  pattern.  The  line 
compass  operator  is  used  over  the  surface  moving  from  left  to  right.  The  location  of  the  peak 
value  of  the  resultant  provides  the  hop  start/stop  time.  The  difference  between  the  true  and  the 
estimated  starting  time  is  the  estimation  error.  It  is  evaluated  in  terms  of  points  of  the  time  axis. 
For  each  SNR  20,  realizations  are  used  with  a  128-point  hop  interval. 

Figures  5.8a  -  5.8c  plots  the  mean  square  estimation  error  (MS)  as  a  function  of  SNR  and 
Daubechies  wavelet  length.  At  an  SNR  2:  10  dB,  the  MS  is  200  or  better.  We  notice  that  the 
shorter  wavelet  (i.e.  Daub-2)  has  the  better  performance  (i.e.,  at  an  SNR  ^  6  dB  the  MSB  is  about 
200).  This  observation  agrees  with  the  notion  that  longer  wavelets  relative  to  shorter  ones  have 
better  frequency  localization  but  have  poorer  time  localization. 
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Figure  5.8a  Hop-Timing  Mean  Square  Estimation  Error  for  Daub-2. 
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6.  Wavelet-Based  Hopping  Time  Detection 


6.1  INTRODUCTION 

As  discussed  earlier,  spread  spectrum  communications  schemes  have  received  ever 
increasing  attention  over  the  past  two  decades  as  numerous  civihan  apphcations  have  joined  mihtaiy 
applications  [22,27,28,32,50,64].  Two  main  assumptions  typically  foimd  in  the  literature  are  that 
the  hop  timing  is  constant  and  known,  and  that  the  hopping  frecjuencies  are  selected  from  a  known 
class  of  candidate  frequencies.  Even  when  the  hop  timing  is  not  assumed  known,  it  is  still  usually 
assumed  constant  [28].  These  assumptions  generally  restrict  the  detection  and  estimation  schemes 
to  frequency  hopping  (FH),  one  of  the  more  popular  spread-spectrum  communications  techniques. 

The  primaiy  goal  of  this  section  is  to  provide  a  new  approach  for  the  detection  and  estimation 
of  frequency  hopping  signals  which  makes  none  of  the  restrictive  assumptions  listed  above.  By  not 
making  such  restrictive  assumptions,  it  is  hoped  that  a  secondary  goal  is  obtained.  This  goal  is  the 
application  to  the  detection  and  estunation  of  other  spread-spectrum  communications  techniques, 
such  as  direct  sequencing,  time  hopping  and  hybrids  of  the  three. 

Section  6.2  briefly  reviews  the  defuoition  of  the  temporal  correlation  function  used  as  the 
backbone  of  the  analyzing  scheme.  Section  6.3  introduces  the  preprocessing  tools  used  to  increase 
the  robustness  of  the  analyzing  to  noise.  Next,  Section  6.4  briefly  explains  how  wavelet  analysis  fits 
in  our  procedure.  Section  6.5  presents  the  detection  scheme  developed.  Section  6.6  presents  the 
overall  detection  algorithm  and  simulation  results.  Finally,  Section  6.7  provides  conclusions  and 
proposed  extensions. 
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6.2  TEMPORAL  CORRELATION  FUNCTION 


The  temporal  correlation  function  (TCF),  also  called  ICF,  of  a  signal  x(k)  is  defined  as: 


TCFj,k,x)=xik^x)x\k-z), 


where  k  is  the  absolute  center  time  and  r  is  the  lag  time,  e^ressed  in  number  of  samples.  Consider 
the  following  analytical  firequency  hopping  signal  x(k)  given  by: 


for  0<k<T,  where  is  the  time  of  the  hop  (or  change  in  frequency)  firom  f  to  f^,  and  where  u(k) 
is  the  imit  step  function.  The  resulting  TCF  function  is  defined  as: 


=rcF,(ifc,-c)+rcFj(it,T)+rcF„(jfc,T), 


where  TCF,(k,z),  TCFjOCjT),  and  TCFjjCk,^  represent  the  1*  2”^,  and  3"*  non  overlapping  terms 
contained  in  the  TCF  expression.  Note  that  computing  the  TCF  of  the  real  fi'equency  hopping  signal 
has  drawbacks  as  additional  “crossterms”  are  present  in  the  resulting  expression,  making  the 
fi-equency  identification  process  more  complex  [63].  Thus,  we  only  consider  the  analytical 
fi-equency  hopping  signal.  In  practical  situations,  the  analytical  signal  can  easily  be  genCTated  by 
applying  a  Hilbert  transform  to  the  real  signal  [50]. 
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Figure  6.1a  presents  the  phase  plot  of  an  analytical  frequency  hopping  signal  x(k)  for  some 
arbitrary y)  and_^,  hoping  time  Tj,Qp=208,  and  for  positive  t  values.  The  combinations  of  the  different 
shifted  versions  of  the  unit  step  functions  force  the  TCF  to  take  on  non-zero  values  only  within  the 
overall  triangular  within  the  regions  shown  in  Figure  6. la.  Note  that: 

1)  TCF,(k,r)  is  a  function  off,  and  ronly.  The  second  term,  TCF2(k,z),  is  a  function  off 
and  ronly,  while  the  last  term,  TCFj2(k,r),  is  a  function  off,f,  k,  and  r. 

2)  The  frequency  hopping  time  T^op  is  located  where  the  region  covered  by  TCFj(k,z)  ends 
and  the  region  covered  by  TCF2(k,T)  begins. 

3)  For  a  given  value  of  r,  the  terms  within  the  triangular  regions  (i.e.,  the  regions  where 
TCF,(k,r)  and  TCF2(k,z)  are  defined)  are  constant,  although  at  different  levels,  while  the 
phase  behavior  within  the  “cross-terms”  region,  TCF,2(k,z)  is  linear.  This  fact  is  further 
illustrated  in  Figure  6.1b  which  plots  the  unwrapped  phase  of  the  TCF  function  for  the  value 
?^30.  It  is  important  to  realize,  however,  that  the  phase  values  over  these  three  regions  are 
a  function  off  znAf  and,  therefore,  not  predictable  without  knowing y)  and^,  which  in 
general  we  do  not.  Nevertheless,  for  any  given  value  of  the  lag,  r<  7^,  this  region  of  cross- 
terms  is  centered  on  foe  hop  time,  T,,^,  another  fact  which  is  exploited. 

6.3  PREPROCESSING  STEPS 

The  main  idea  behind  foe  proposed  scheme  is  to  take  advantage  of  foe  TCF  phase  behavior 
along  foe  time  axis  k  (i.e.,  for  fixed  values  of  As  shown  earlier  foe  unwrapped  TCF  phase  along 
foe  time  axis  is  constant  prior  and  after  foe  frequency  jump,  while  it  is  linear  in  a  region  centered 
around  foe  frequency  hopping  time,  resulting  in  a  constant-ramp-constant  phase  behavior  along  foe 
time  axis  k,  as  illustrated  in  Figure  6.1b.  Differentiating  such  phase  leads  to  a  pulse  centered  around 
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the  hopping  time,  as  shown  in  Figure  6.1c.  Detecting  the  edges  of  the  pulse  is  then  all  what  is 
needed  to  identify  the  hopping  time,  as  it  can  then  be  estimated  as  the  midpoint  between  the  two 
edges.  Therefore,  the  hopping  time  detection  problem  can  be  viewed  as  an  edge  detection  problem, 
which  the  wavelet  transform  is  well-matched  to  address. 

However  the  additive  noise  contained  in  the  communication  signal  results  in  phase  noise, 
degrading  the  quality  of  the  resulting  pulse.  Furthermore,  the  differentiating  step  increases  the  effects 
due  to  noise.  Thus,  we  apply  median  filtering  before  and  after  the  differentiating  step  to  minimize 
these  phase  noise  effects.  The  main  advantage  behind  this  filter  is  that  it  preserves  the  ramp  and  step 
behavior  and  eliminates  outliers.  The  length  of  the  median  filter  operation  was  selected  in  order  to 
preserve  the  step  discontinuities  present  in  TCF(k,  r)  for  fixed  r.  Further  details  may  be  found  in 
Overdyk  [63]. 

6.4  WAVELET  TRANSFORM  BASED  DETECTION 

Edge  detection  is  an  important  problem  in  numerous  applications  ranging  firom  image 
processing  to  transient  detection,  and  wavelet  transforms  have  been  used  extensively  for  detecting 
discontinuities  in  a  given  signal  or  its  derivatives.  Recall  that  wavelets  may  be  used  to  detect 
discontinuities  in  a  signal  or  its  derivatives,  if  the  chosen  wavelet  function  is  able  to  represent  tiie 
highest  order  derivative  present  in  the  signal  function,  as  any  wavelet  with,  at  least  p  vanishing 
moments,  can  be  used  to  detect  a  discontinuity  in  the  (p-1)**  derivative  [46].  For  our  problem,  we 
are  interested  in  identifying  pulse  edges,  i.e.,  a  signal  discontinuity,  and  the  Haar  wavelet  is 
sufficient  for  the  task.  In  addition,  detecting  discontinuities  can  be  done  quite  simply  using  the 
decimated  version  of  the  wavelet  transform  [46]. 

The  presence  of  noise  makes  identification  of  discontinuities  more  complicated.  In  such  a 
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C3.S6,  the  sverEging  of  several  scales  can  enhance  the  wavelet’s  ability  to  detect  discontinuities  in 
noise.  The  basic  idea  is  that  for  true  discontinuities,  the  spikes  will  line  up,  while  those  associated 
to  noise  will  not  [46].  As  a  result,  we  averaged  the  first  two  scales  in  our  simulation  to  improve  the 
robustness  of  the  detection  scheme  to  noise  degradations. 

Figure  6.2a  plots  the  result  obtained  by  averaging  the  first  two  scales  of  the  DWT  applied 
to  the  unwrapped  phase  of  the  TCF  expression  for  an  analytical  fi-equency  hopping  (FH)  signal.  The 
FH  function  has  frequencies  on  either  side  of  the  hop,  located  at  time  sample  208, =  6.250  MHZ 
and_^  —  22.727  MHZ.  The  SNR  level  is  equal  to  10  dB.  Figure  6.2b  illustrates  the  resulting  wavelet 
transform  of  the  pre  processed  phase  of  the  TCF  obtained  at  lag  ?=30.  Figure  6.2b  shows  that  the 
wavelet  transform  clearly  detect  the  location  of  the  pulse  ends,  as  expected.  Note  that  the  width 
between  each  spike  obtained  firom  taking  the  wavelet  transform  of  the  TCF  phase  along  the  time  axis 
k  (i.e.,  for  fixed  lag  time  r)  increases  as  r  increases.  The  next  step  sums  all  the  values  which 
represent  the  edges  of  the  cross-terms  in  the  TCF  in  45°  and  135  “directions  so  that  they  reinforce 
each  other,  as  illustrated  in  Figure  6.3b.  Figure  6.3c  shows  the  “detection  vector”  obtained  fi^om  this 
summation.  Further  details  may  be  found  in  Overdyk  {63].  Note  that  the  resulting  spike  is  located 
at  the  hopping  time  Th^=2C8. 

6.5  DETECTION  SCHEME 

Once  the  detection  vector  has  been  formed,  a  decision  must  be  made  as  to  whedier  or  not  a 
hop  has  occurred  within  the  fimne.  E3q>eriments  suggested  that  the  variance  of  the  detection  vector 
would  be  the  best  indicator  of  whether  or  not  a  hop  had  occurred.  As  a  result,  the  threshold, 
is  chosen  as  a  multiple,  k,  of  the  variance  of  the  detection  vector  whai  no  hop  has  occurred  within 
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the  frame.  The  threshold  determination  was  also  guided  by  the  fact  that  the  cost  associated  with  the 
probability  of  a  missed  detection,  P„  =  [1  -  probability  of  detection  (Pd)]  far  outweighs  the  cost 
associated  with  the  probability  of  a  false  alarm,  Pf^,  as  the  hopping  time  estimation  is  only  the  first 
step  in  a  complete  frequency  hopping  signal  detection  scheme.  This  is  described  in  detail  in  the  next 
paragraph. 

Note  that  once  the  hopping  times  are  estimated,  the  signal  frequencies  need  to  be  extracted 
to  demodulate  the  actual  message.  This  step  can  easily  be  done  by  appl3dng  frequency  analysis  to 
the  estimated  hopping  intervals.  Thus,  in  the  case  of  false  alarms,  frequency  analysis  would  show 
the  same  frequencies  in  two,  or  more,  consecutive  hopping  intervals,  resulting  in  no  message 
degradation.  However,  a  missed  hopping  time  will  result  in  degradations  in  the  frequency  estimation 
step,  and  errors  in  the  decoded  message.  Receiver  operating  characteristics  (ROC)  curves  were 
generated  for  six  SNR  levels  and  an  appropriate  threshold  chosen  for  each.  Further  details  regarding 
the  choice  of  specific  thresholds  may  be  found  in  Overdyk  [63]. 

6.6  DETECTION  ALGORITHM  AND  RESULTS 

As  stated  earlier,  it  was  desired  to  make  as  few  assumptions  as  possible  on  the  nature  of  the 
frequency  hopping  signal.  With  this  goal  mind  the  assumptions  were  limited  to  three: 

1.  Known  spread  spectrum  frequency  range.  This  range  was  assumed  to  be  limited  from 
1  MHZ  to  24  MHZ  in  the  simulations  conducted, 

2.  Known  minimum  hopping  time,  T^^ This  parameter  was  chosen  to  be  equal  to  256 
sample  points  in  our  simulations.  At  a  sampling  rate  of  50  MHZ,  this  translates  into  a  minimum 
hopping  time  of  5.12  ps, 

3.  Minimum  frequency  dififerential,  Af,  for  frequency  hops.  This  parameter  was  chosen  to 
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be  1  Khz  in  our  simulations. 


Detection  Algorithm 

Using  the  tools  described  above,  the  algorithm  steps  for  the  detection  and  estimation  of 
frequency  hopping  signals  in  noise  can  now  be  enumerated  as  follows: 

1.  Transform  the  real  signal  into  an  analytic  signal. 

2.  Segment  data  into  frames  of  length  less  than  or  equal  to  the  minirmiTn  hopping  time, 

’^hopjnin-  assumption  ensures  that,  at  most,  one  hop  will  be  present  in  the  processing 
frame. 

3.  Compute  the  temporal  correlation  Junction  on  each  frame. 

4.  Extract  the  TCP  phase  information  and  unwrap  the  phase  along  the  time  axis  k. 

5.  Apply  a  median  filter  to  the  phase  of  the  TCP  along  the  time  axis  k,  of  length  5.  This  step 
is  done  to  reduce  the  noise  effects  prior  to  differentiating,  since  differentiating  accentuates 
the  effects  due  to  noise. 

6.  Differentiate  the  phase  information  along  the  time  axis  k.  This  step  changes  the 
unwrapped  phase  of  the  TCP  from  a  ramp  function  to  a  pseudo-pulse. 

7.  Apply  a  second  median  filter  along  the  time  axis  k  of  lenjgth  25.  The  length  of  25  has 
proven  to  work  well  with  the  chosen  for  our  simulations  [63].  This  step  is  done  to 
again  remove  the  effects  due  to  noise  which  were  accentuated  by  the  differentiation  opaation 
in  step  6. 

8.  Calculate  DWTs  along  the  time  axis  k  (i.e.,  for  each  lag,  t,  of  the  TCP)  using  the  Haar 
wavelet.  This  step  is  done  to  detect  the  discontinuities  at  the  edges  of  the  cross-terms  region. 
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9.  Average  the  wavelet  coefficients  of  the  &st  two  scales  of  the  DWT. 

10.  Perform  a  45°/135°  summation  across  all  values  of  lag,  t,  to  obtain  a  detection  vector 
which  has  time  as  its  index. 

1 1.  Threshold  the  data  in  the  resulting  detection  vector  obtained  in  step  10.  When  detected 
above  the  flireshold,  the  maximum  peak  value  time  index  represents  the  estimated  hopping 
time. 

Simulation  results. 

Simulations  were  conducted  to  test  the  effectiveness  of  the  detection  and  estimation 
algorithm  given  above.  Five  himdred  trial  experiments  were  conducted  in  six  different  signal  to 
noise  ratios  (SNR)  between  -3  to  15  dB.  The  basic  idea  behind  the  experiments  was  to  simulate 
signals  that  had  already  been  segmented,  as  specified  in  steps  1  and  2  of  the  algorithm.  The  problem 
then  becomes  to  determine:  a)  whether  or  not  a  firequency  hop  exists  within  the  given  firame;  b)  to 
estimate  the  hopping  time  when  a  firequency  hop  is  detected. 

Communication  signals  were  generated  by  choosing  random  hopping  time,  T^^p,  and  hopping 
frequencies  selected  randomly  in  a  predefined  range.  The  resulting  signal  is  a  signal  with,  at 
most,  one  hop  which  can  be  from  any  frequency,/),  to  any  hopping  frequency,  such  that  1  MHZ  ^ 
/;/2^24MHZ. 

Table  6.1  shows  the  probability  of  detection,  P^,  the  probability  of  false  alarm,  Pp^,  and  the 
percentage  of  errors  m  classification  for  the  selected  threshold,  T^ahou*  ^  of  the  six  SNRs 
considered.  Note  the  entries  under  the  column  labeled  “A”  represents  a  multiple  of  the  variance  of 
the  detection  vector  generated  from  a  “no  hop”  fimne.for  each  respective  SNR  level.  The  colmnn 
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labeled  “%  Error”  shows  the  percentage  of  mis  classifications  (i.e.,  the  percentage  of  false  alarms 
plus  misses).  Note,  also,  that  the  low  probability  of  false  alarm  was  sacrificed  for  higher 
probabilities  of  detection  for  reasons  discussed  earlier.  For  example,  the  entries  in  the  row  for  SNR 
=  3  dB  show  that  if  a  1 1.4  %  mis  classification  rate  and  a  Pfa=0-1961  can  be  tolerated,  then  we  can 
expect  to  detect  89.53%  of  the  hops  in  a  given  frequency  hopping  signal 

Simulation  results  are  shown  in  Table  6.2.  The  column  labeled  “Avg.  Error”  gives  the 
average  error  obtained  at  each  SNR  level.  For  example,  at  the  SNR  level  of  3  dB,  out  of  all  the  hops 
which  were  detected,  the  average  distance  firom  the  true  hopping  time  was  10.48  sample  points  (i.e., 
4.1%  of  the  minimum  hopping  time).  Columns  with  numeric  headings  indicate  the  hop  detection 
probability  within  a  given  percentage  of  For  example,  at  the  SNR  level  of  3  dB,  the  column 

labeled  “1%”  indicates  that  36%  of  all  detectedhoT^s  were  located  within  1%  of  7;^^  ^  or  within  2 
points  of  the  true  hop  time,  Similarly,  72%  of  all  detected  hops  were  located  within  5%  of 
Thopjnin  or  within  12  points  of  the  true  hop  time, 

6.7  CONCLUSIONS 

This  section  considered  the  application  of  temporal  correlation  functions  and  wavelet  analysis 
to  the  detection  and  estimation  of  fi-equency  hopping  signals  in  additive  white  Gaussian  noise.  The 
algorithm  developed  has  only  two  restrictive  assumptions:  1)  a  minimum  hopping  time;  2)  a 
minimum  fi-equency  differential.  Thus,  it  can  find  applications  where  tiie  minimum  hopping  time 
is  not  held  constant;  i.e.,  in  time  hopping  signals  and  hybrid  techniques  involving  txihsc  frequency 
hopping  or  time  hopping.  We  showed  how  the  detection  problem  can  be  formulated  as  an  edge 
detection  problem,  and  how  wavelet  analysis  can  be  used  in  the  hopping  time  detection  problem. 
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Third,  we  introduced  preprocessing  techniques  designed  to  improve  the  robustness  of  the  detection 
and  estimation  scheme  in  noisy  environments.  Results  show  that  the  scheme  is  robust  to  noise 
degradations  down  to  3  dB. 

Although  the  algoritiim  developed  succeeded  in  meeting  its  goals,  improvements  can 
potentially  be  obtained  by  taking  advantages  of  the  specific  two-dimensional  structure  of  the 
hopping  pattern.  Thus,  a  possible  extension  involves  considering  the  problem  as  an  image 
processing  or  pattern  recognition  problem,  due  to  the  specific  triangular  pattern  generated  by  the 
TCP  of  frequency  hopping  signals.  As  a  result,  applying  more  sophisticated  two-dimensional  edge 
detection  schemes  may  improve  the  robustness  of  the  detection  and  estimation  scheme  to  noise 
distortions.  As  in  the  first  part  (i.e.,  chapters  1-5),  improvements  in  the  estimation  of  the  TCP  (or 
ICF)  will  improve  the  performance  of  the  estimation  procedure,  that  is  improve  detection  and 
estimation  of  hop  times.  Such  extensions  are  left  for  further  research. 
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Figure  6.1;  The  unwrapped  phase  information  of  ( a  )  the  TCF,  (  b  )  for  T!=30,  and  ( c  )  after  i«fP«»t-i.nf{artng 
and  median  filtering  for  x=30. 
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Figure  6.2:  DWT  coefiBdents  computed  for  eachvalue  of  lagx(a)onthepreprocessedphaseof  the 
TCF(t,t),  and  for  t=30. 
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Figure  6.3:  Detection  vector  constructed  by  performing  a  457135°  summation.  ( a  )  shows  effects  of  the 
summation  in  (  b  )  where  flie  arrows  indicate  the  direction  of  summation.  ( c  )  shows  a  plot  of  the  actual 
detection  vector. 


SNR 

Pfe 

%  Error 

15  dB 

140 

1 

0 

0.0 

10  dB 

30 

0.9866 

0.0196 

1.4 

6dB 

15 

0.9844 

0.1569 

3.0 

3dB 

11 

0.8953 

0.1961 

11.4 

OdB 

0.3529 

20.4 

SEIDHHI 

10.6927 

0.3333 

31.0 

Table  6.1:  Detection  statistics  of 500  experiments  applying  the  detection  and  estimation  algorithm. 


SNR 

(dB) 

1% 

5% 

10% 

15% 

20% 

30% 

40% 

50% 

75% 

100% 

Avg.  Error 
(#  of  samples) 

15 

0.790 

0.984 

0.992 

0.996 

0.996 

0.998 

1.00 

1.00 

1.00 

1.00 

2.22 

10 

0.726 

0.964 

0.974 

0.978 

0.982 

0.986 

0.986 

2.70 

6 

0.558 

0.926 

0.940 

0.950 

0.970 

0.970 

0.970 

5.46 

3 

0.360 

0.720 

0.758 

0.794 

0.828 

0.862 

0.874 

0.882 

0.886 

0.886 

10.48 

0 

0.116 

0.302 

0.418 

0.510 

0.572 

0.684 

0.752 

0.768 

0.796 

0.796 

28.48 

-3 

0.174 

0.276 

0.382 

0.456 

0.568 

0.614 

0.646 

0.686 

30.99 

Table  6.2:  Estimation  statistics  for  the  500  experiments  at  each  SNR  level  using  the  detection  and  estimation 
algorithm  described  in  Section  B.  Columns  with  numeric  headings,  show  the  probability  that  estimated  hops 
are  found  within  a  given  distance,  expressed  in  percentage  of  Thopjmn^  from  true  hopping  times. 
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7.  CONCLUSIONS  AND  RECOMMENDATIONS 
7.1  CONCLUSIONS 


Wavelet  analysis  of  2-D  correlation  functions  is  a  new  area  of  investigation.  It  can  be 
applied  to  the  interception  of  communication  signals.  This  work  aims  at  applying  wavelet 
analysis  to  the  instantaneous  correlation  function  to  identify  frequency  hopped  signals.  The 
instantaneous  correlation  function  (ICF)  of  the  complex-valued  FH  signal  is  shown  to  have  a 
cellular  (diamond)  pattern,  where  each  hop  generates  one  main  diamond  structure.  Inside  a  given 
diamond,  the  ICF  of  the  signal  consists  a  single  complex  exponential  component  representing  the 
hop  frequency  in  the  delay  direction  and  some  noise.  In  the  time  direction,  inside  a  given 
diamond  pattern,  the  wavelet  transformed  data  tends  to  be  a  constant  perturbed  by  noise.  The 
intersections  of  the  diamond  with  the  time  axis  determine  the  hop  start/stop  times  while  the 
width  of  a  given  diamond  corre^onds  to  the  hop  interval. 

The  wavelet  transform  of  the  ICF  surf^e  generates  a  number  of  surfaces.  The  wavelet 
surface,  at  any  scale,  emphasizes  the  frequency  hops  which  reside  there  and  attenuates  spectral 
components  that  do  not  belong  to  the  particular  scale  (i.e.,  bandpass  filter)  under  consideration. 

If  we  apply  the  wavelet  transform  along  the  delay  axis,  we  can  address  the  interception 
problem  in  two  different  ways.  We  can  visually  inspect  the  wavelet  surfaces  to  identify  and 
classify  based  on  the  structure  the  modulation  due  to  an  FH  signal.  This  also  allows  to  obtain  a 
rough  estimate  of  the  hop  time  interval.  Alternatively,  we  can  apply  a  processing  scheme  which 
can  be  used  to  automate  the  interception  task.  This  processing  scheme  estimates  the  hop 
start/stop  times,  the  hop-scale  pattern,  and  the  hop  frequency. 

The  estimation  of  the  hop  start/stop  times  can  be  addressed  using  an  edge  detection 
approach.  We  applied  a  compass  operator  to  find  edges  in  the  wavelet  filtered  ICS.  The  hop  scale 
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pattern  is  obtained  by  applying  an  energy  analysis. 


The  frequency  of  each  hop  can  be  extracted  from  ttie  wavelet  surface  at  the  proper  scale, 
or  from  the  original  time  data  using  the  hop  time  parameters. 

Visual  inspection  of  the  wavelet  smfaces  permits  the  identification  of  FH  signals  at  SNR 
levels  of  3  dB  and  above.  Other  modulation  schemes  such  as  ASK,  PSK,  and  MFSK  will  only 
have  cellular  patterns  on  one  of  the  wavelet  surfaces,  that  is  their  frequency  bandwidth  does  not 
span  more  than  one  wavelet  scale.  Hop  timing  estimation  shows  that  the  hop  start/stop  times 
can  be  estimated  with  an  accuracy  of  12  to  17.5  per  cent  at  SNR  levels  of  6  dB  or  better. 

The  performance  of  longer  duration  wavelets  is  better  than  that  of  shorter  ones  since 
longer  wavelets  have  better  spectral  energy  concentration  than  shorter  ones.  The  success  rate  of 
frequency  estimation  from  the  wavelet  surfaces  showed  that  the  probability  of  correct  frequency 
estimation  from  the  wavelet  surface  is  1.0  for  input  SNR’s  of  0  dB  and  above.  A  frequency 
estimation  success  rate  of  1.0  requires  an  SNR  level  of  -3  dB  or  about  -3  dB  to  -5  dB  using  the 
time  data  directly  or  wavelet  surfaces,  respectively.  The  minimum  SNR  required  for  an 
automated  estimation  of  hop  times  is  6  dB. 

Processing  along  frie  time  axis  allows  detection  of  the  transition  times  of  frequency 
and  time  hopped  signals.  The  current  implementation  is  limited  to  work  with  one 
transition  per  observation  interval  but  permits  robust  detection  at  an  SNR  of  3  dB. 


7.2  RECOMMENDATIONS  FOR  FUTURE  WORK 

There  are  other  ways  to  define  or  estimate  the  instantaneous  (or  temporal)  correlation 
fimction.  A  study  to  evaluate  different  candidate  ICS  (or  TCS)  should  be  performed. 
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For  processing  of  the  ICS  along  the  delay  axis  we  recommend  the  following: 

i)  Refine  the  automatic  recogmtion  of  the  cellular  structure  of  the  FH  signal  over  the  wavelet 
surfaces. 

ii)  Improve  the  performance  of  the  hop-scale  identification  at  lower  SNR  by  reexamining  the 
equali2ation  of  die  spectrum  of  the  ICF  surfaces. 

iii)  Investigate  other  wavelet  types,  and  the  use  of  other  definitions  for  the  instantaneous 
correlation  fimction. 

iv)  Combine  information  fi*om  different  wavelet  surfaces  to  improve  parameter  estimation. 

For  processing  the  ICS  along  the  time  axis  we  recommend  investigation  of  extensions  for 
successful  operation  when  more  than  one  transition  (jump)  is  present  during  the  observation 
interval. 

Finally,  both  ^proaches,  transformation  over  time  and  delay,  will  benefit  firom  an 
improved  edge  detection  scheme.  This  problem  should  be  addressed  in  more  detail. 
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